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We apply the unbiased weak-coupling continuous time quantum Monte Carlo (CTQMC) method 
to review the physics of a single magnetic impurity coupled to s-wave superconducting leads de- 
scribed by the BCS reduced Hamiltonian. As a function of the superconducting gap A, we study 
the signature of the first order transition between the singlet and doublet (local moment) states 
on various quantities. In particular we concentrate on the Josephson current with to tt phase 
shift, the crossing of the Andreev bound states in the single particle spectral function, as well as 
the local dynamical spin structure factor. Within DMFT, this impurity problem provides a link 
to the periodic Anderson model with superconducting conduction electrons (BCS-PAM). The first 
order transition observed in the impurity model is reproduced in the BCS-PAM and is signalized 
by the crossing of the low energy excitations in the local density of states. The momentum resolved 
single particle spectral function in the singlet state reveals the coherent, Bloch-like, superposition 
of Andreev bound states. In the doublet or local moment phase the single particle spectral function 
is characterized by incoherent quasiparticle excitations. 
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I. INTRODUCTION 

Magnetic degrees of freedom in superconducting envi- 
ronments have attracted considerable interest due to the 
underlying competing effects. Already a classical spin 
oriented along the z-axi4^ embedded in a superconduct- 
ing host generates a localized state within the supercon- 
ducting gap. As a function of the interaction strength 
this excitation crosses the Fermi energy thereby trigger- 
ing a first order transition between a ground state with 
vanishing total electronic spin and a ground state with 
nonzero total electronic spin. 

For a quantum spin, the Kondo effect sets in. Being 
a Fermi surface instability, the opening of the supercon- 
ducting gap competes with Kondo screening and ulti- 
mately leads to a local moment regime. This transition 
is accompanied by a to tt phase shift in the Josephson 
current. In the local moment regime the 7r-shift occurs 
since a Cooper pair tunneling through the junction nec- 
essarily accumulates a phase tjUS^. 

The interest in the problem has been renewed in the 
last decade by the rapid progress in nanotechnology 
which made a direct experimental realization of quan- 
tum dots coupled to superconducting leads feasible so 
that many experiments have been designed to directly 
measure the to tt transition of the Josephson current. 
Experiments using a carbon nanotube^ ^ ^ but also InAs 
nanowireJ^ as a quantum dot coupled to superconduct- 
ing leads were able to observe the sign change of the 
Josephson current by increasing the gate voltage and 
thus manipulating the number of electrons on the quan- 
tum dot. The effect of the changing electron number 
the behavior of such systems has been extensively 
studiecP^^^^^ and the theoretical expectation of the 
collapse of the Kondo effect if the superconducting gap A 



exceeds the Kondo temperature Tk has been confirmed 
by experiments of Buitelaar et al.^^. 

From the numerical point of view, a combination of 
algorithmic development and computational power has 
allowed for a more detailed study of the proble m us- 
ing the numerical renormalization grou quan- 
tum Monte Carlo simulations^^ ^ ^^ as well as functional 
renormalization group calculation^^. Most numerical 
works present in the literature only present either the 
study of the Josephson current^^ or the study of 

the spectral properties of the Quantum dot^^. One of the 
goals of this article is to use the weak coupling CTQMC 
method^^ to compute the Josephson current as well as the 
spectral functions for the same parameter set in order to 
present a comprehensive study of the to tt transition of 
a Josephson quantum dot. Our numerically exact data 
clearly confirms the picture of a first order phase transi- 
tion from a singlet phase linked to the 0-j unction regime 
of the Josephson current to a doublet phase correspond- 
ing to the TT-j unction regime. 

In addition to numerical efforts, many analytical ap- 
proximations have been introduced to tackle different as- 
pects of the physics of the problem. The non crossing ap- 
proximation has been used to show that Andreev bound 
states crossing the Fermi energy are connected to the 
to TT transition of the Josephson current^^. Perturbative 
methods as well as mean field theory have brought a quite 
complete understanding of the phase diagram featuring 
the and tt phases as well as the intermediate phases 0' 
and tt'^^ Another method employed by several au- 

thors is the introduction of different analytically solvable 
effective models, which are valid in different limit J^^'^^'^. 
These models are very useful to acquire an intuitive un- 
derstanding of the physics. We will present the study of 
an effective Hamiltonian for the limit of a superconduct- 
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ing gap A much larger than the bandwidth to support 
the interpretation of the CTQMC data. 

Another motivation of this paper, is to study within 
dynamical mean field theory (DMFT)^^ the periodic 
Anderson model with an s-wave BCS-conduction band 
(BCS-PAM). Within this approximation, the BCS- 
PAM maps onto the single impurity Anderson model 
with superconducting baths supplemented with a self- 
consistency condition. We will show that the physics 
of the impurity model can be taken over to the lattice 
case. In particular the first order transition observed in 
the impurity model is reproduced in the BCS-PAM and is 
signalized by the crossing of the low energy excitations in 
the local density of states. The momentum resolved sin- 
gle particle spectral function in the singlet phase reveals 
the coherent, Bloch-like, superposition of Andreev bound 
states. In the doublet or local moment phase the single 
particle spectral function is characterized by incoherent 
quasiparticle excitations. We provide an understanding 
of this in terms of models of disorder. 

The paper is organized as follows. After introducing 
the model in Sec. [TTj we discuss in Sec. |III| an effective 
toy model valid in the limit of a superconducting gap, 
A, much larger than the bandwidth W. This simple toy 
model goes a good way at understanding certain aspects 
of the underlying physics. A brief outline of the employed 
CTQMC result including the proof of Wick's theorem for 
each configuration in the Monte Carlo simulation will be 
presented in Sec. [IV| The results of the toy model are 



then compared to the results of the CTQMC simulation, 
which are discussed in detail in Sec. [Vl Sec. IVII is dedi- 
cated to the study of the BCS-PAM within DMFT. We 
include an appendix |A] featuring the proof of a general 
determinant identity needed for the proof of Wick's the- 
orem for every configuration in the CTQMC. 



II. MODEL 

The physics of a quantum dot coupled to two supercon- 
ducting leads (L=left, R=right) via a hybridization term 
is captured by the single impurity Anderson model with 
the leads described by the BCS mean-field Hamiltonian: 



(1) 
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The operators ^ ^ are creation operators for electrons 
with a z-component of the spin a and momentum k in 
lead a, dj. is a creation operator of an electron with a 
2;-component of the spin a on the quantum dot. ^/^ = 
e{k) — /i = —2tcos{k) — /i is the dispersion relation for 
the electrons in the leads, where we assume, that the 
dispersion is independent of the lead index a, and = 
Cd — fi is the position of the dot level. Throughout this 
paper, we will express all quantities in units oft = 1. The 
superconducting order parameter has a modulus A and 
a phase The parameter V characterizes the strength 
of the hybridization, and U corresponds to the Coulomb 
blockade. 

Since the Hamiltonian does not conserve the electron 
number as a consequence of the BCS-term, we use the 
standard trick of rewriting the Hamiltonian in terms 
of creation and annihilation operators of quasiparticles, 
which for spin up are identical to the electrons, but cor- 
respond to holes in the spin down sector. This can also 
be expressed as a canonical transformation: 



d 



(3) 

Using the new operators, the Hamiltonian can be writ- 
ten in a Nambu notation: 
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k,a 



(4) 



with Hu = —U{d^^d^ — — the Nambu spinors 



d = ( 1^ I , Ck,a 



^k,i,a 



(5) 



the matrices 



Ea(fc) 
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and the Pauli matrix 

'1 



-1 r 



(7) 



For practical reasons, we use the following definition 
for the single particle Green's function throughout Sec. 
HI] to Sec. El 

13 



(8) 



With this definition, the resolvent operator G'^{iuJm) = 
(— zcjml — Hq^) can be used to obtain the Green's 
function of the noninteracting system: 



{-iuj, 

Yl 
N 



E 



(9) 
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III. EFFECTIVE HAMILTONIAN IN THE 
LIMIT A/iy ^ oo 

To gain a deeper understanding of the physics on the 
quantum dot, it is useful to search for analyticahy solv- 
able toy models. We will study an effective model, which 
reproduces the physics of the Hamiltonian ([T]) in the limit 
A/VF oo, where W is the band width. To derive the ef- 
fective model, we look at the limit A ^ oo of the Green's 
function in Eq. ([9| . The superconducting order parame- 
ter A appears only in the matrix Eq;(/c), thus we examine 
the behavior of this matrix for large values of A. This 
can easily be done by diagonalizing Eq,(A;) for (j)^ = 0* 



Running of eigenenergies 
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Ua. (10) 



Let us first look at the limit A ^ oo of the transformation 
matrix Ua^ which for brevity is not a unitary matrix. 



A 



1 1 
-1 1 



(11) 



The diagonal matrix in Eq. (10) can be considered in a 



similar manner and we obtain for lim Ec^(/c) = Eq 

A^oo 



-A 
A 



-A 
-A 



(12) 



Using this result, for large values of A the sum over k 
and a in Eq. ([9| can be carried out yielding 

(13) 

This is exactly the free Green's function obtained from a 
Hamiltonian of the form: 

iJeff = -\/2y(cV^d + dV^c) + c^EooC + dUdd + Hu- 

(14) 

i^eff describes a system consisting of one bath site c con- 
nected by a hybridization term to the correlated quantum 
dot d. The dispersion of the bath has completely van- 
ished, as the superconducting band gap becomes much 
larger than the bandwidth. 

We chose a basis of the 16 dimensional Hilbert space 
and write the Hamiltonian as a matrix, which subse- 
quently can be diagonalized. As we have restricted 
the parameter space for the Monte Carlo simulations to 
e^^ = and /i = in the original Hamiltonian of Eq. Q, 
we will use the same parameters for the exact diagonal- 
ization results. 



Ground state of the effective model 



The ground state of the system (14) can be deter- 
mined by diagonalizing the Hamiltonian i^eff- As de- 
picted in Fig. [l] the energy levels cross at a critical 



FIG. 1: (Color online) Eigenenergies of the effective Hamil- 
tonian (14) for varying U. The fixed parameters are given by 
y = 0.5 and A = 1. The crossing of the two lowest levels is 
clearly seen at U ~ 1.7. The ground state for [/ < 1.7 is a 
singlet state. For larger values of [/, the twofold degenerate 
doublet state becomes energetically more favorable. 



Running of ground state weights 




FIG. 2: (Color online) For A < 1.412 the ground state is the 
singlet state from Eq. ( p^ . If A is increased, the weight a of 
the single occupied states | i ) and | i, T ) decreases in favor of 
the states with a double occupied quantum dot, corresponding 
to the weights [3 and 7. At A = 1.412 the ground state 
changes to the twofold degenerate doublet state given in ( 16 ) 



and the weight of the states with a single occupied quantum 
dot b increases with A. The parameters in this plot are V — 
0.5 and U = 1.0. 



value oiU = Uc and a similar behavior can be observed 
by varying A with a corresponding critical value Ac- 
For U < Uc and A < Ac, the ground state is given 

by \^s) = -ce(in,o)-|o,n))-/3(iT,i) + u,T))- 

7(U,i) + I T,T)), with the notation 4 I 0,0) = |cr,0) 
and cZj. |0,0) = |0,cr). Note, that we are using the un- 
physical basis introduced in Eq. (|3|. To interpret this 
ground state it is better to return to the physical ba- 
sis by inverting the canonical transformation in Eq. ([3| 
and transforming the vacuum state | 0, 0) | i, i )• The 
ground state can then be rewritten in the physical basis 



4 
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Double occupancy 




FIG. 3: (Color online) Double occupancy ( J | J | J | J | ) of the 
quantum dot in the effective model at /3 = 200 and V = 
0.5. This plot can be understood as a phase diagram of the 
effective model, as the phase boundary is accompanied by a 
sharp decay of the double occupancy. 



very sharp drop of the double occupancy on the phase 
boundary can be observed, which evolves to a jump at 
T = 0. Here the larger values of the double occupancy 
are connected to the singlet phase, while the lower values 
belong to the doublet phase, where single occupancy is 
favored. This can be understood by studying the expec- 
tation value of the double occupancy in the ground state. 
In the singlet phase, we obtain 



and for the doublet phase: 



(17) 



(18) 



From the behavior of the weights 7 and a shown in 
Fig. [2] it is clear that the double occupancy increases 
with A in the singlet phase and decreases in the doublet 
phase. 

Note, that many of the results presented in this pa- 
per can be observed either at fixed U or A as can be 
conjectured from Fig. [3j 



\^s)=a{\U)-\ll)) 

+ /3(|0,fl) + |fl,0)) +7(10,6) 



in,n: 



(15) 



This state is clearly a singlet state, corresponding to 
a Kondo singlet between the quantum dot and the bath 
with the dominant weight a. The states representing 
a pairing on the quantum dot or in the bath have the 
suppressed weights P and 7 for small values of A but 
grow more important if A is increased as is shown in 

Fig.m 

At U > Uc^ the ground state changes and we 
get the twofold degenerate ground states | t ) ~ 

a (I n, T ) - 1 n, i )) (I T, n ) + 1 1, n )) and 1 ^Pd,i ) = 

a(|0,T)-|0,J,)) + 5(U,0) + I T,0)), rewritten in the 
physical basis: 

I Vd.T ) = a (1 f , ) - 1 T, n )) + 6 (1 0, f ) + 1 n, f )) 
I V'rf.i ) = a (1 1, 6) - 1 1, n )) + 5 (1 0, 1 ) + 1 n, I )) . 

(16) 

This two-fold degenerate ground state has a z-component 
of the total spin ±1/2 and hence corresponds to a local 
moment. 

B. Phase diagram 



C. Proximity effect 

To gain further insight in the sign change of the local 
pair correlations ( J| J| f ^ ' ^"'^H we calculate the ground 
state expectation value of the local pair correlations in 
the effective model (14). For the singlet phase, we obtain 



\d]d\ 



i's) 



(/?lfI,fI)+7lfL0)) 
2Re(/3*7) > 0. 



(19) 



Clearly, only terms describing the pairing on the quan- 
tum dot contribute to the pair correlations, whereas the 
Kondo singlet of electrons on the quantum dot and in the 
bath does not. From Fig. [2j it is obvious that the result- 
ing pairing correlation is positive and increases with A. 
This illustrates the proximity effect, as a pair field in the 
bath induces a pair field on the quantum dot. 

On the other hand, in the doublet phase, we obtain 



\d\d\ 



<0. (20) 



As in the singlet phase, only the states corresponding 
to a pairing on the quantum dot contribute to the pair 
correlations. The local moment part of the ground state 
does not generate pair correlations. As the weight a in 
the doublet phase ground state is positive and decreases 
with A (see Fig. |2|, the local pair correlations have a 
negative sign in contrast to the positive sign in the singlet 
phase and decrease with A. 



To further illustrate the phase transition between the 
singlet state | ips ) and the doublet states | il^d,u )' dou- 
ble occupancy ) of the quantum dot in the ef- 
fective model is shown in Fig. [3] At low temperature a 



D. Spectral function 

Using the Lehmann representation, the spectral func- 
tion A||(a;) of the effective model is easily calculated. It 
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Spectral function 



Dynamical spin structure 
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FIG. 4: (Color online) Spectral function A^^{lj) of the effec- 
tive model for different values of A at /3 = 200, [/ = 1 and 
V = 0.5. The ^-peaks have been broadened by a Gaussian 
function of width a — 0.04 for better visibility. 
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FIG. 5: (Color online) Dynamical spin structure factor S{(jj) 
of the effective model at p — 200. The phase transition from 
the singlet-phase to the doublet-phase for U — 1 and V — 
0.5 occurs at A 1.412. At this point a transition from 
a gapped excitation to a peak at a; = corresponding to 
a local magnetic moment in the doublet phase is observed. 
To visualize the ^-functions, a Gaussian broadening of width 
a — 0.05 has been applied. 



is defined by 



71,171 



( n I J I I m ) 



(21) 
The 



with the matrix elements Mnm 
spectral function is shown in Fig. [4j Comparing this plot 
to the numerical solution of the full model as depicted in 
Fig. [l4j we observe, that the simple model already shows 
the important feature of an excitation at the position 

= at the critical value of A. Even though for very 
small values of A, the Kondo resonance at c<; = can 
not be seen in the simple model, we see a precursor of 
the Kondo resonance as a pole of the Green's function, 
which develops into a resonance if we increase the number 
of sites in the batlP^. 

A careful analysis reveals, that the low frequency sig- 
nature of the spectral function reflects the excitation be- 
tween the two lowest lying states of the spectrum. These 
states are the ground states of the singlet and the dou- 
blet phase and therefore, the position uj of the excitation 
marks precisely the energy difference of the two ground 
states. At the critical value of A = 1.412, the level cross- 
ing occurs and leads to a vanishing energy difference of 
the two ground states, meaning that the excitation be- 
tween the two states lies now precisely at = 0. 



E. Dynamical spin structure 

Like the spectral function, the dynamical spin struc- 
ture factor S{u;) can be calculated using the Lehmann 



representation: 



(n|5+|m) 6{uo^En-E^). (22) 



In the Monte Carlo simulation, a numerically more sta- 
ble quantity is obtained by replacing 5+ by Sz in the 
above equation. This quantity is completely equivalent 
to S{uj), as we only make use of the 5'[/(2)-symmetry of 
the problem, and is therefore used in the following. 

In the representation (22) of 5'(a;), it is clear that the 
dynamical spin structure factor will show excitations at 
frequencies corresponding to the energy needed to flip the 
spin on the quantum dot. Therefore, the dynamical spin 
structure factor is very well suited to determine whether 
the system is in the singlet or in the doublet regime. 

In Fig. |5]the phase transition from the singlet phase to 
the doublet phase is reflected by the fact, that in the sin- 
glet phase, a gapped excitation can be observed, whereas 
in the doublet phase, a peak at = emerges, which 
corresponds to a local magnetic moment. 



F. Dynamical charge structure 

The dynamical charge structure factor N{u;) can be 
defined by the Lehman representation 

iVH = -| ^ |(n| n - <5„,™ |m)|' e-'^^™5(w+i?„-i?™). 

n,m 

(23) 

As for the other spectral functions, the charge struc- 
ture factor N{ij) shown in Fig. [6j exhibits a sharp change 
of its behavior at the phase transition for the critical 
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Dynamical charge structure 



N{uj) 




FIG. 6: (Color online) Dynamical charge structure factor 
N{oo) of the effective model at /3 = 200. We have used the 
same parameters as for Fig. [4] 



value of the superconducting gap A. We observe, that 
the charge structure shows a finite gap for all values of 
A and that for large values of A, the gap increases in a 
slightly nonlinear manner. 

A more detailed study of the matrix elements con- 
tributing to the charge structure factor reveals, that be- 
cause of correlation we have completely different exci- 
tations than for the spectral function. In fact, the most 
prominent excitations are excitations from the respective 
ground states in the two different phases to higher energy 
states with structure similar to that of the ground states. 



IV. CTQMC 
A. Basic outline of the algorithm 

For the numerically exact solution of the BCS- 
Anderson-model, we used the weak coupling CTQMC- 
methocP^, which is based on a perturbation expansion 
around the limit of U = 0. Following the presentation of 
the CTQMC-algorithm in^^, we will shortly outline the 
basic principles of the method. 

As pointed out iipHMl the interacting Hamiltonian Hu 
in Eq. Q can up to a constant be rewritten as 



Hrj = 



U 



=±1 



a^)(d\d^-al) (24) 



introducing the parameters to minimize the sign prob- 
lem. For the present case, a choice of = a^^ = ^ -\- sS 
with 6 = 1+0+ was found to completely eliminate the 
sign problem at half filling, even after the complex phase 
factors exp(i^Q,) in the Hamiltonian were introduced. 
Using perturbation theory, the partition function Z of 



the full Hamiltonian Q can be written as: 



Zo 



X (T (nT(Ti) - af ) . . . (ni(r„) - a^) )o • 

with the number operators ha = d^^da and the thermal 
expectation value ( • )o = ^ Tr [e~^^°»]. As Hq is a 
noninter acting Hamiltonian, %ick's theorem holds, and 
the expectation value {T{h^{ri) — a|) . . . {hi{rn) — a^))^ 
can be cast in a determinant of a matrix Mc^ of size 
2nx2n, where Cn is a configuration of vertices {r^, Si}. In 
contrast to the formulation for the Hubbard model given 
in^^, we do not need to include an index for the lattice 
site as we only have one correlated site, the impurity. The 
Matrix Mc^ is not block diagonal for the two spin sectors 
in the case A 7^ 0, so we cannot factor the determinant in 
two determinants of n x n matrices. Finally, the partition 
function of the model is given by 



Z 



E 



U 



detMc„, 



(26) 



where the sum runs over all possible configurations Cn 
of vertices as irP^. The matrix Mc^ is defined by 



M 



/GO,(ri,ri)-ai 



V GO,(ri,r„) 



(27) 



using the 2x2 Green's function matrices G2^(r, r') 

iTd\ir)d,ir'))„ iTd\ir)d,ir'))J ^"^^ ^^^^ " \^ 

A similar reasoning yields an expression for the thermal 
expectation value (0(r)) = ^ Tr [e-^^O(r)] of the full 
model: 



(o(t); 



E 



. 2 J 



^detMc„((0(T)))c„ 



Ec„(f) detMc„ 



(28) 



Here {{0{r)))cr, is the contribution of the configuration 
Cn to the observable 0(r), which is given by 



((0(r)))c„ 



(r(nT(Ti)-ai)...(nx(T„)-apO(T))o 
(T(nT(Ti)-a|)...(nx(T„)-ap)o 



(29) 

Both, the numerator and the denominator of the above 
Eq. ( 29 ) can be written as determinants of matrices using 
Wick's theorem. Eq. (28) is the central relation of the 



CTQMC algorithm, because starting from this equation, 
the Metropolis-Hastings- Algorithm can be employed to 
generate a Markov chain of configurations Cn- At this 
point, we have to interpret (^) detMc^ as the statis- 
tical weight of a given configuration Cn what in general 
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is impossible, as detMc^ is a complex number. There- 
fore, we have to replace (^) detMc^ by its modulus 
and account for the phase in the measurement of the 
observables. Fortunately, in the present case, the statis- 
tical weights are always real and nonnegative, so that we 
can simply calculate the contribution to the observable 
0(r) for a given configuration Cn in the Markov chain 
as {(0(r)))c„. 



tity of interest, the calculation of the contribution 
((T7I71/ . . . 7j^7^/))c'^ is tedious and time consuming. 
Luckily for every configuration Cn a relation similar to 
Wick's theorem can be found, which greatly simplifies 
the calculation of higher Green's functions. It is closely 
connected to the determinant identity ( A2 ) proven in ap- 



pendix |Aj The application of the ordinary Wick's theo- 
rem to the denominator and the numerator of Eq. (29) 
yields 



B. Wick's theorem for each configuration 

For the measurement of higher Green's functions 
of the form (T7j7i/ . . . 7j^7^/ ), where 7/ stands for 
c^i,(^i,meas) or 4, (^^,meas) depending ou the quan- 



((^7l7l' •••7m7m'))c^ 



det Bc^ 
detMcJ 



(30) 



where we have defined the matrix Bc^ G 

(Q(2n+m) X (2n+m) 



( 



M, 



(T7ld(ri))o . 



(T7ld(T„))o . 
(Tdt(ri)7i')o ... (Tdt(r„)7i')o k^^hv)^ ■ 



V{rdt(ri)7„.'>o 



r 



{r7t,d(r„))o 



(31) 



Defining the matrices B^ G (C^'^^^^y^^^^\ we can 
make use of the determinant identity ( A2 ) 



(T7jd(r„))o 
V{Tdt(ri)7i'>o...(rdt(r„)7.')o {Tlh')o J 



yielding 

detBc^ _ 1 

det Mc^ ~ (detMc^)^ 



(32) 



/detB^i 



det 



\det B^^ 



detB^^^ 



det B^^ , 



(33) 

From Eq. 29 it is obvious, that det B^ / det M^^ is iden- 
tical to the contribution of the configuration Cn to the 
one particle Green's function {T^j^if). Hence, Wick's 
theorem holds for every configuration Cn and is given by 



n{Tlhi'))c„ ... ((T7l7i'))c„ 
det : ■.. : 



(34) 



This relation is particularly useful in a simulation mea- 
suring multiple physical observables as measurements of 



single particle Green's functions can be reused in an eco- 
nomic way. 

V. NUMERICAL RESULTS 

In this section, we present the results obtained by 
CTQMC simulations for the model ([T]). We restrict our- 
selves to the case of half filling, e^^ = and /i = 0. In the 
first part of this section, we will discuss the results for 
static quantities including the Josephson current, double 
occupancy and pair correlations on the quantum dot. We 
then proceed to dynamical quantities such as the single 
particle spectral function and the dynamical spin struc- 
ture factor. 



A. Josephson current 

The Josephson current fiowing through the Quan- 
tum dot can be calculated directly within the CTQMC 
method, as it is given by an equal time Green's function: 



{ja ) = i-X^ Yl - dlck,a,a ) (35) 

k,a 

We show here our results for the Josephson current at 
an inverse temperature of /3 = 50 as a function of the su- 
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Josephson current for difTerent values of A 



Josephson current for different values of A 




FIG. 7: (Color online) Josephson current in the junction 
regime 






FIG. 9: (Color online) Josephson current in the tt junction 
regime. 



perconducting gap A. For small values of A, we observe 
a sinusoidal form of the Josephson current as a function 
of the phase difference (j) with increasing amplitude, as A 
increases (see Fig. [t]). 



Temperature dependence of the Josephson current 



Josephson current for different values of A 
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FIG. 10: (Color online) Josephson current for different tem- 
peratures. The current phase relations do not intersect at one 
single point as suggested by the NRG results of Karrasch et 

aim. 



FIG. 8: (Color online) Josephson current in the 0^ and tt' 
junction regime. 

This parameter regime is known as the 0- Junction 
regime, because the Josephson current /j((/>) = §^ has 
a zero with positive slope at ^ = 0, corresponding to a 
minimum in the grand potential Q at (j) = (see Fig. 5 
in referenc^^. 

If the value of A is further increased, the behav- 
ior of the Josephson current changes, as in the region 
A 0.15 . . . 0.35 the Josephson current shows a zero be- 
tween ^ = and (j) = TT. (see Fig. [8|. This leads to a 
minimum in the grand potential at tt and the parame- 
ter regime is called 0' or tt' regime depending on which 
minimum of the grand potential is the global on^^. The 
behavior of the Josephson current is in accordance with 
the behavior of the double occupancy seen in Fig. [l2j 
as in the same parameter region, where we observe the 
0' to tt' transition, the drop of the double occupancy as 
a function of (j) can be observed, which is linked to the 
change of the curvature of the current-phase relation of 
the Josephson current. 



For larger values of A, the sign of the Josephson cur- 
rent changes and the grand potential shows now a single 
minimum at = tt, this regime is therefore called the tt 
regime, (see Fig. [9|. 

The picture for the behavior of the grand potential 
as a function of (j) that we get from the current phase 
relation of the Josephson current agrees very nicely with 
the results presented by Benjamin et al.^. 

The current phase relations for the different phases pre- 
sented here were also extensively studied by Karrasch et 
al. using the fRG and NRG methods^^, Choi et al. us- 
ing the NRG methocPS, as well as by Sia no and Egger 
using; the Hirsch-Fye QMC method^SEHHI^ g^^^ though 
the numerical exactness of certain results has been de- 
bated, the results of all numerical works show very good 
qualitative agreement and are confirmed by the present 
results. 

In the literatur^^l^, the temperature dependence of 
the current phase relation of the Josephson current has 
been discussed. We show CTQMC results in Fig. 
which look very similar to the Siano and Egger result^. 
As CTQMC is numerically exact, our result suggests that 
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the crossing of all curves in one single point^ at 
found in the approximate finite temperature NRG is not 
universal. 



7, =0 



B. Double occupancy 



Double occupancy of the Quantum dot 



We learned from the toy model described in Sec. [IIT 
that the system exhibits a phase transition from the sin- 
glet phase to the doublet phase as U is increased. This 
picture is consistent with the NRG results of Bauer et 
dl^. The phase transition can be observed in the double 
occupancy {n^fii ) of the quantum dot, which is propor- 
tional to where O is the grand potential. At T = 0, a 
sharp step function of the double occupancy is expected. 
While the T = regime is not directly accessible to quan- 
tum Monte Carlo calculations, we calculated the double 
occupancy for different temperatures using the CTQMC- 
method. The results are shown in Fig. [TT] From the 
data, it is obvious that with decreasing temperature the 
curves converge to the step function of the limit T = 0, 
which is a clear sign for a first order phase transition, 
reflecting a level crossing of the two ground states. This 
is in complete accordance with the results for the toy 
model. 



Double occupancy 
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FIG. 11: (Color online) Double occupancy (h^hi) of the 
quantum dot at A = 1.0. The data shows a jump in the 
double occupancy becoming sharper with decreasing temper- 
ature. 

It is interesting to correlate the Josephson current as a 
function of the phas e diff erence (j) = (j)L — (j)R ioi various 
values of A (see Sec. V A), with the double occupancy on 



the dot. As depicted in Fig. 12, for very small values of 
A as well as for A >^ 0.4, we see that the double occu- 
pancy is a constant function of (j). This corresponds to a 
current-phase-relation for the Josephson current fixed in 
either the tt- or the 0-j unction regime. For intermediate 
values of A, we observe a far more interesting behavior of 
the double occupancy: At a certain value of ^, the dou- 
ble occupancy drops to a smaller value. This drop is of 
course smeared out by the finite temperature, but can be 






FIC. 12: (Color online) Double occupancy of the quantum 
dot as a function of the phase difference = — for 
different values of A. 



understood as a way to drive the phase transition from 
the 0- to the 7r-j unction regime by the phase difference 0. 



C. Pair correlation 

In agreement with the NRG result of Choi et al.l^ as 
well as with the mean field results by Salkola et. al.l^, 
we obtain the local pair correlation on the quantum dot 
shown in Fig. 13 For small A, the local pair correlation 



increases because of the proximity effect, as an increas- 
ing magnitude of the pair field A in the leads induces a 
growing pair correlation on the quantum dot. The sharp 
sign change at the critical value of A observed at zero 
temperature is smeared out at finite temperatures, but 
the qualitative behavior is exactly the same as for the 
effective model discussed in Sec. IIII CI We therefore con- 
clude, that the sign change of the pair correlation is due 
to residual pairing on the quantum dot in the doublet 
phase which decreases with A. 

The same qualitative behavior of the local pair corre- 
lation is al so obs erved, if U is changed instead of A as 
discussed in^^. The sign change of the local pair cor- 
relation A^^ is traditionally expressed as a 7r-phase shift 
in Ad. 



D. Spectral function 

All quantities studied so far suggest that a first order 
phase transition occurs when we tune the system from 
the 0- Junction to the tt- Junction regime. This can be 
confirmed by studying dynamical quantities such as the 
spectral function. 

In Fig. 14 we show the spectral function A{lj) of the 
quantum dot as a function of A. The data has been cal- 
culated from the CTQMC data for the Green's function 
G^^(r) using stochastic analytic continuatioiP^I^. This 
method works especially well for the low energy spectrum 
and sharp excitations while the high energy spectrum and 
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Local pairing correlation 



Spectral function 




FIG. 13: (Color online) Local pair correlation = {^\^\) 
as a function of A. We observe the same behavior as Choi 
et al!^, which is also in very good agreement with the pair 
correlation expected for the effective model discussed in |III C| 



excitation continua are more difficult to resolve. Inside 
the gap, the formation of Andreev bound states can be 
seen very well. 




FIG. 14: (Color online) Spectral function A{(jj) as a function 
of A for the parameters P = 100, [/ = 1.0 and V = 0.5 at 
half filling and zero phase difference between the two super- 
conductors. 



E. Dynamical spin structure factor 



In the region of A ^ we see the Kondo-resonance. As 
a function of growing values of A and as a consequence of 
the opening of the quasiparticle gap at the Fermi level, 
the Kondo resonance evolves to Andreev bound state. 
Note that at the mean-field level, the Kondo resonance 
merely corresponds to a virtual bound state. Opening a 
quasiparticle gap at the Fermi level drives the lifetime of 
the this virtual bound state to infinity. In the parameter 
region which corresponds to the 0- Junction regime of the 
Josephson current (A ^ 0...0.1), we observe Andreev 
bound states with excitation energies approaching co = 0. 
This corresponds to the crossing point in Fig. 14 and 



has also been observed by Bauer et al. for fixed A and 
increasing U in^. 

The comparison of the Quantum Monte Carlo data 



shown in Fig. 14 with the result obtained from the ef- 
fective model discussed in Sec. |IIID| is particularly in- 
sightful. The spectral signature is very similar except for 
the lack of the Kondo resonance due to the finite size 
of the effective model. In the effective model, the An- 
dreev bound state excitation corresponds to the energy 
difference between the ground states of the singlet and 
the doublet phase. The position A at which the Andreev 
bound states cross at cc; = has been identified as a clear 
sign for the crossing of the ground states of the singlet 
and doublet phases. Hence, we interpret the crossing of 
the Andreev bound states in the CTQMC data as a very 
strong sign for a level crossing and hence a first order 
phase transition from the singlet to the doublet phase in 
the full model. 



In addition to the spectral function, the dynamical spin 
structure factor S{uj) defined in Eq. 22 , provides a way 
of characterizing the phases of the system. For A = 0, we 
clearly see a suppressed spectral weight at = and a 
peak which corresponds to the characteristic energy scale 
of the Kondo temperature Tk- From the peak position, 
we obtain a rough estimate for the Kondo temperature 
of Tk ^ 0.06. 

From A ^ 0.05 onwards, spectral weight is accumu- 
lated at a; = ultimately forming a pronounced sharp 
local moment peak for large values of A. As the Kondo 
temperature is a measure for the energy required to break 
the Kondo singlet, we expect the Kondo effect to break 
down at a value of A ^ Tk- This is indeed observed in 



Fig. 15 



The signature of the breakdown of the Kondo reso- 
nance also shows up in the spectral function plotted Fig. 



14 Since the Kondo resonance stems from a screening 
of the magnetic moment by conduction electrons in an 
energy window Tk around the Fermi level, the opening 
of a single particle gap of order Tk destroys the Kondo 
resonance giving way to an Andreev bound state. 

The breakdown of the Kondo resonance is accompa- 
nied by a change of the curvature in the current-phase- 
relation of the Josephson current which is a precursor for 
the transition to the 0' phase (see the curves for A = 0.05 
and A = 0.08 in Fig. u\. We also observe that after the 
transition from the tt- to the tt- regime has occurred 
(see the current-phase-relation of the Josephson current 
of Fig. [8| the peak at finite co vanishes and all the spec- 
tral weight is accumulated in the very sharp local moment 
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peak at = 0. 



Chargegap as a function of A 



Dynamical spin structure factor 




FIG. 15: (Color online) Dynamical spin structure factor S{(jj) 
as a function of A for the parameters p — 100, U — 1.0 and 
V — 0.5 at half filling and zero phase difference between the 
two superconductors. For A = we can roughly estimate 
the Kondo-Temperature Tk ~ 0.06 from the peak position of 
5H. 



F. Charge gap 

From the dynamical charge structure factor, we can 
determine the gap Ac to local charge fluctuations on the 
dot with two different methods^^ . One way to extract the 
charge gap is to read off the peak position of the lowest 
lying excitation in the dynamical charge structure factor 
obtained from the charge correlation function Cc(r) = 
(n(r)n) — (n) (n) via stochastic analytic continuation. 
The other way of extracting the charge gap from Cc(r) is 
based on the fact, that the charge structure factor N{uo) 
is linked to Cc(r) via 



(36) 



If N{(jo) is sharply peaked around a certain value cjp, we 
can approximate N{(jo) by N{(jo) ^ d{uj — ujp). This cor- 
responds to Cc{r) ^ e~'^^p. Therefore, a least squares 
fit of an exponential function e~^^p to Cc{t) in a region 
where only one single mode dominates, can reveal the 
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FIG. 16: (Color online) Charge gap Ac as a function of A. 
We calculated the dynamical charge structure factor from the 
charge-charge correlation function Cc{r) using stochastic an- 
alytic continuation and extracted the charge gap using two 
different methods. First, we read off the charge gap directly 
from the stochastic analytic continuation data, secondly, we 
calculated the charge gap from the charge-charge correlation 
function. The straight line is a linear fit through the numeri- 
cal data. 



frequency cOp at which N{u;) is peaked. The applicability 
of the method can be seen in the half logarithmic plot 
of Cc(r), where a sharply peaked charge structure factor 
N{u;) is reflected by a region, in which Cc{t) can be well 
approximated by a straight line. 

The data obtained using these methods is shown in Fig. 
[161 In the context of the effective model discussed in Sec. 



IIIF[ we observe, that the behavior of the charge gap of 
the full model clearly differs from that of the effective 
model. Especially, we do not see any signature of the 
phase transition in the behavior of the charge gap. 

The charge gap opens approximately linearly with A. 
It is very hard to extract the charge gap from the numer- 
ical data at small A, therefore we can only extrapolate 
to A = 0. Here, it appears, that we have a finite charge 
gap even in the absence of superconductivity. 

The fact that the local charge fluctuations remain 
gaped confirms the picture that the to tt transition 
occurs only in the spin sector. 



VI. DMFT 

A. Periodic Anderson Model with BCS conduction 

band 

In the previous sections, we have studied the first or- 
der phase transition in the impurity model ([T]). As the 
dynamical mean field theory (DMFT) provides a link be- 
tween impurity models and lattice models, we can ask 
the question if the singlet to doublet phase transition 
observed in the impurity model is also realized in a cor- 
responding lattice model. 
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An appropriate lattice model will of course include a 
U{1) symmetry breaking term like the impurity model 
([T]) does, and in fact in the framework of the DMFT, 
a periodic Anderson model extended by the BCS mean 
field Hamiltonian (BCS-PAM) for the conduction band 
electrons corresponds to the impurity model presented 
in the previous sections^. The Hamiltonian of the BCS- 
PAM is given by: 



H = He + Hf + Hy 



(37) 



with 



Hc = J2 C(^)4,.5fc.. - A E {iAki + h-c.) (38) 
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(40) 



k,a 



We have considered a square lattice with hopping matrix 
element t between the conduction electrons such that: 



^(/c) = —2t {cos{kax) + cos{kay)) . 



(41) 



Note, that the impurity model ([T]) has a large range 
of applications in the DMFT ranging from the attractive 
Hubbard model with U{1) symmetry broken solutions 
studied in references'^ to the BCS-PAM, which is con- 
sidered here. 

The treatment of this model within DMFT involves the 
same steps as for the impurity model ([T]), introducing a 
particle- hole transformation for the spin down operators. 
The Hamiltonian can then be cast in the form H = Hq + 
Hu with 



k k k 

(42) 

and Hu = (^i/,T ~ I) {^ifd ~ k)' 

have used the same Nambu-spinor notation as in Sec. [n| 
with the exception, that d operators have been renamed 
/ to be consistent with the literature^^ 



B. DMFT with superconducting medium 

The standard DMFT can be easily adapted to a su- 
perconducting bath using the Nambu formalism^^. We 
obtain the self consistency equation for a finite lattice 
with N sites expressed by a 2 x 2 matrix equation: 



G*f(iw„) 



k 



(43) 



FIG. 17: Double occupancy of the / sites in the BCS-PAM. 
In the proximity of the critical value of [/, we observe two 
different solutions of the DMFT self consistency cycle. The 
upper (red) branch is generated, if we start the DMFT algo- 
rithm with a self energy E = 0, while we obtain the solution 
shown by the lower (blue) branch if we take the self energy of 
the data point at = 0.44 as the starting point of the DMFT 
iterations. 



(Tf(r)ft) is the fun 



(3 

Here, G^{iuJn) = —Jdre~ 



Matsubara Green's function of the reference model, 
(icOn) is the Matsubara /-Green function of the bare 



lattice model and 5] is the self energy. Equation (43) 



can be solved by iteration starting usually at a self en- 
ergy 5]^ = 0. From G^{iuJn)^ the bare Green's function 
Goi'^^n) of the reference model, can be calculated using 
Dyson's equation Qq = G^ + X)^. The reference 
model, which is now described by Qq and the interaction 
part of the Hamiltonian can subsequently be solved us- 
ing the CTQMC method yielding G^^iuOn) for the next 
DMFT iteration. 



C. Hysteresis 

In the DMFT, we can calculate the double occupancy 
(/| ih,if\ if 14 ) ^f /-sites, which is together with the 
assumption of a homogeneous system proportional to ^ . 
Therefore, we expect a jump in the double occupancy to 
appear at a critical value of /7, if we have a first order 
phase transition as in the impurity problem. 



Figure 17 shows our result for the double occupancy of 
the / sites as a function of U. Depending on the initial 
choice of the self energy in the DMFT cycle, we obtain 
two different solutions. If we start with the local Green's 
function of the bare lattice model, which corresponds to 
a self energy E = 0, we obtain the upper branch of the 
hysteresis. The lower branch is obtained by taking the 
self energy of the solution in the strong coupling phase 
at /7 = 0.44 as starting point for the DMFT cycle. The 
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coexistence of two solutions is a strong hint that a first 
order phase transition occurs. 

It should be noted that beginning at a value of U ^ 
0.34, the upper branch of the hysteresis becomes unsta- 
ble, i.e. the inherent fluctuations of the Monte Carlo 
results suffice to drop from the upper branch of the hys- 
teresis to the lower branch after a certain number of iter- 
ations. Increasing the number of Monte Carlo measure- 
ments delays the drop to the lower branch to a higher 
number of iterations. This behavior can be understood 
in the following way: In the coexistence region, the grand 
potential Q of the upper and lower branch of the hystere- 
sis cross at a certain value of U. For small values ofU, Q 
is minimal on the upper branch, while the lower branch 
is metastable, for larger values of /7, however, the stable 
solution is the lower branch. 

In the strong coupling phase and on the lower branch 
of the hysteresis, the Monte Carlo results suddenly de- 
velop a finite magnetization corresponding to a frozen 
spin. This is due to divergent autocorrelation times in 
the Monte Carlo simulation and is linked to the physical 
formation of a local moment. 



D. Local dynamical spin structure factor 

To further classify the weak and strong coupling 
phases, we calculate the local dynamical spin structure 
factor S{ij) = ;^ 'S'(q, cj). The Lehmann representa- 
tion for S{u;) is given by Eq. (22), where in this case 



weak coupling regime and shows a characteristic energy 
scale required for flipping a spin. 

The lower branch of the hysteresis represents the 
strong coupling phase and shows a clear local moment 
peak in the dynamical spin structure factor at a; = 0. 

This behavior reflects exactly the single impurity 
physics discussed in the previous section where we ob- 
served the Kondo effect in the weak coupling phase and 
the formation of a local moment in the strong coupling 
phase. 



E. f-Density of states 

In order to investigate the behavior of the /-bands at 
the phase boundary and to be able to compare with the 
single impurity model, we calculate the density of states 
for the /-sites pff directly from the local Green's function 
G(r) using the stochastic analytic continuation method 
for different values of U. From Fig. [T9j one can recog- 
nize the signature of the impurity physics (see Sec. VD), 
namely the crossing of Andreev bound states in the vicin- 
ity of the first order transition at /7 ^ 0.35. Note, that 
we have only shown the level crossing for the impurity 
model if A is changed, but for varying [/, the crossing of 
the Andreev bound states in the impurity model ([T]) has 
been observed by Bauer et al.^^. Clearly in the lattice 
model, one expects the Andreev bound states to acquire 
a dispersion relation which shows up as a finite width in 
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FIG. 19: Density of states for the f-electrons as a function 
of U for the parameters V — 0.5, A = 2, /i = e/ = and 
13 - 100. 



FIG. 18: Dynamical spin structure factor for the upper and 
the lower branch of the hysteresis in Fig. [T7| Clearly, the up- 
per branch of the hysteresis corresponds to a singlet solution, 
while the lower branch shows a local moment. 



F. Dispersion relation of Andreev bound states 



As in the impurity case, S{(jo) is a measure for the en- 
ergy needed to flip the spin on an /-site. Figure [18] shows 
the result for the local dynamical spin structure factor on 
both branches of the hysteresis. The solution correspond- 
ing to the upper branch of the hysteresis is linked to the 



We have seen in the previous subsections, that the lo- 
cal physics of the single impurity model can be carried 
over to the lattice case within the DMFT approximation. 
Here, we concentrate on unique features of the lattice 
model (37), namely the dispersion relation of the f-bands 
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as obtained by analyzing the single particle spectral func- 
tion. 

Using the local self-energy of the DMFT, this 
quantity is extracted from the Green functions 



(44) 



and 



(45) 

where G^J^''(zcJn), G^^f (zcj^), ^^^^ij^^n), ^^^ij^^n) de- 
note the noninteracting Green functions for the corre- 
sponding orbitals in the unit cell. 

Using the stochastic analytic continuation, these 
Green's functions can be rotated to real frequencies, 
yielding in principle the spectral function A(k, cj). For 
each k-point and real frequency this quantity is a 4 x 4 
matrix since we have a 2 x 2 Nambu spectral function 
for each combination of / and c orbitals. Our analysis of 
the spectral function is based on the basis independent 
quantity A{k,uj) = Tr A(k, cj). 



Trace of the spectral function A(k,a;) at C/ = 0.125 
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FIG. 20: Trace of the spectral function A{\<i,uj) at /3 = 100 
in the singlet regime. The parameters of the simulation were 
given by [/ = 0.125, V 0.5, A = 2 and /x = e/ = 0. 



Fig. 20 plots this quantity in the singlet phase. The 
overall structure of the spectral function is similar to the 
structure observed for the bare BCS-PAM characterized 
by the four bands: 



±(k) = ±^JV^ + ^2(k)/2 ± E(k)VV2+£:2(k)/4 

(46) 

where £^(k) = ^/eH^)~^r~^ . The bands with dominant 
c-character, £^±(k) = £^^'+(k), at high frequencies are 
well separated from the bands of dominant /-character 
at low frequencies, £^£(k) = E±-{k). For the considered 
bare parameters, V is the smallest scale and sets the 
magnitude of the dispersion relation of the /-band. In 
particular expanding in V gives: 



Em = ± 



o 



(- 



i?(k) V^(k) 



(47) 



Starting from the point of view of the impurity model, 
which as seen above accounts very well for overall form 
of the k-integrated /-spectral function, £^£(k) may be 
perceived as the dispersion relation of the Andreev bound 
states. 



Trace of the spectral function A(k, a;) at [/ = 
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Trace of the spectral function A(k, a;) ai U — 0.125 
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Trace of the spectral function A(k, u) dX U = 0.2 
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Trace of the spectral function A(k, u) aX U = 0.275 
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FIG. 21: Trace of the spectral function A(k,cj) at /3 = 100 
in the singlet regime for increasing interaction U . The width 
of the /-bands clearly decreases and the dispersion becomes 
weaker. The parameters of the simulations were given by 
V = 0.5, A = 2 and /X = e/ = 0. 



The singlet phase is continuously connected to the 
/7 = point. Starting from this limit, we can account 
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for the Hubbard U within a slave boson approximation'^ 
which win renormahze the hybridization matrix element 
to lower values. Owing to Eq. [47|this suppresses the dis- 
persion relation of the /-electrons. This aspect is clearly 



observed in Fig. 21 



Trace of the spectral function A(k,a;) at U — 0.5 
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Trace of the spectral function A(k, u) at U = 0.55 
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FIG. 22: Trace of the spectral function A{k,uj) at /3 = 100 
in the doublet regime for different values of U. Here, we only 
show the /-bands. The parameters of the simulation were 
given hy V = 0.5, A = 2 and /i = e/ =0. 



now reads: 

(3 (3 

Smf = JdrJ dr'i\r)g-\r-ry{ry 





UrUz 



drf^(r)f(7 



(49) 

where ft = (/|,/|) and g{ T — t') corresponds to the 
bath Green function. To account for disorder, the z- 
component of the f-spin is sampled from the box distri- 
bution rriz G [—Mz^Mz]. Averaging over disorder at each 
iteration in the DMFT cycle yields the spectral function 
shown in Fig. [23j As apparent, the disorder average 
generates a finite lifetime. 



Disorder average of the trace of the spectral function A(k, a;) 
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FIG. 23: Trace of the spectral function A(k,cj) as obtained 
from using Eq. [49] for the impurity action. The z-component 
of the local moment is sampled from the box distribution 
rriz G [—Mz^Mz]. The parameters used for this plot were 
given by y = 0.5, U = 0.5, A = 2 and = 0.0375. Here, 
the calculations are carried out on the real time axis such that 
no analytical continuation is required. 



In the doublet phase, U > Uc^ the paramagnetic slave- 
boson mean-field approach fails. In this state, the /-spin 
is frozen and in the DMFT cycle we have imposed spin 
symmetric baths thereby inhibiting magnetic ordering. 
The QMC data of Fig. 22 points to a very incoherent 



/-spectral function. It is therefore tempting to model 
this state in terms of spin disorder: the spin of the /- 
electrons on each site is static and points in a random 
direction. To provide some support for this picture we 
stay in the dynamical mean field framework but consider 
a mean-field decomposition of the Hubbard term in the 
action of the impurity problem: 



U 



Urriz 



(^/,T-^/,i) (48) 



This mean field approximation, accounts for the local 
moment formation with z-component of spin m^. The 
corresponding mean-field action of the impurity model 



VII. CONCLUSION 

We have shown that the weak-coupling CTQMC al- 
gorithm is an extremely powerful unbiased tool to com- 
pute thermodynamic as well as dynamical quantities of 
impurity models in superconducting environments. The 
method can cope very well with a complex phase of the 
superconducting order parameter thereby allowing for 
the calculation of the Josephson current. Our detailed 
results for the impurity problem confirm the picture of 
a first order phase transition between a single and dou- 
blet state. It is accompanied by a tt phase shift in the 
Josephson current. Being completely unbiased, our ap- 
proach provides the first numerically exact results for this 
model Hamiltonian. 

Within DMFT, the physics of the BCS-PAM is 
mapped onto the single impurity Anderson model sup- 
plemented by a self-consistency loop. We have shown 
that within this approximation, the physics of the impu- 
rity model can be carried over to the lattice. In particular 
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at fixed superconducting order parameter A tiie first or- 
der transition between a singlet and local moment state 
as a function of growing values of U shows up in a hys- 
teresis behavior of the double occupancy. Furthermore, 
the low energy features of the local /-spectral function 
are reminiscent of the Andreev bound states with van- 
ishing excitation energy (i.e. a crossing point) at the 
critical coupling. Within the DMFT approximation, we 
can look into the single particle spectral function. In the 
singlet phase, the low energy features can be interpreted 
in terms of a dispersion relation of Andreev bound states. 
This state is continuously linked to the U = limit. In 
the doublet state or local moment regime, the low energy 
features of the spectral functions are incoherent. We pro- 
pose to understand this in terms of models of disorder. 
In particular in this state, the spin dynamics of the /- 
electron is frozen and since we are considering paramag- 
netic states it points in a random different direction in 
each unit cell. A simple model of disorder following this 
picture accounts very well for the observed incoherent 
spectral function. 



VIII. ACKNOWLEDGMENTS 

We thank Julia Wernsdorfer for interesting discussions 
and for bringing up the subject in her diploma thesis. We 
also wish to thank Volker Meden for fruitful discussion 
and advice. Part of the calculations were carried out at 
the Leibniz Rechenzentrum in Munich on HLRB2. We 
thank this institution for allocation of CPU time. 

DJL also thanks Jutta Ortloff and Manuel Schmidt 
for many valuable discussions as well as Burkhard Ritter 
for critical reading of the manuscript. FFA would like to 
thank the KITP, where part of this work was carried out, 
for hospitality (NFS Grant PHY05-51164). We thank the 
DFG for financial support. 



APPENDIX A: PROOF OF THE DETERMINANT 
IDENTITY 

In this section a general determinant identity is proven, 
which can be used to derive Wick's theorem for contri- 
butions of a configuration Cn to physical observables. 
Let us define the vectors Ui, Vi G and the num- 
bers aij e C. Further, let A G C^><^ be a ma- 
trix of rank m. We define the non-singular matrices 
Mn e c(^+^)x(^+^) and Ay G C^^+^)xi^+'^) by: 



/A Ui ... Un\ 

vi^ ail ... C^ln 



A- - I ^ 



\Vn^ C^nl ... OinnJ 



(Al) 

With these definitions, the following determinant iden- 
tity holds: 



/ det All 



detMn(det A) 



det 



det Ain^ 



ydet Ani . . . det A^ 



(A2) 

The identity can be proven by induction in n. It is 
trivial for n = 1, so we have to start with n = 2, where 
we have to show 



det M2 det An det A22 det A12 det A21 



det A 



det A det A 



det A det A 



(A3) 



For the following calculations, we introduce several vec- 
tors: 



Vi \ 2 1 
),u =v 



(A4) 















1^22 -1/ 




\ / 



m+2 



(A5) 

Let us define the expanded matrix Cgx of a square matrix 
C as the matrix C expanded by one row and one column 
containing a unit vector: 



C 

0^ 1 



(A6) 



As a last definition, we introduce the abbreviation hij = 
Vi-^A~^Uj. Using these notations, we can write the ma- 



trices Aij as 



(A7) 



To calculate the determinant det Aij , we use the matrix 
determinant lemma det(A+uv-^) = (l+v-^A~-'^u) det A, 
yielding 



det Aij 

det Apx 



1+ vg^(Ae. + ul^v^' )-\i^J (l+v^^ A-iu?,). 

The inverse matrix of (Aex + ^ijV"'" ) can be obtained 
from the Sherman-Morrison formula and a tedious calcu- 
lation making use of the special form of the vectors and 
matrices gives the result 

det Aii 



det A 



= Oil 



(A9) 



From this, the right hand side of Eq. (A3) can be easily 
obtained. For the left hand side, we have to perform 
an analogous calculation using the decomposition of the 
matrix M2: 



M2 



^llex 



11-^ 2 2-^ 



(AlO) 
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Again, we apply the matrix determinant lemma two (A2) holding for n: 
times and insert the Sherman-Morrison formula to calcu- 
late the inverse matrix of (Ang^ + uJ^vJ^ ). Simplifying 
the result as far as possible, we finally arrive at 

det M2 
det A 



= (an - 611) (^22 - ^22)-(q^12 - ^12) (q^2i - ^21) . 



detMn+i(detA)(^-^) = det 



<^2,n+l 



If we compare (All) with (A9), it is clear, that Eq. ( |A3| ) 
holds. 

We now assume that for a certain value n G N Eq. 



(A2) holds. For n + 1, we can cast the matrix Mn+i in a 



form, where we can make use of Eq. (A2) holding for n: 



( A 



Mn+l = 



V 



U2 

<^2,2 

<^n,2 



\^n+l <^n+l,2 . . 



<^2,n+l 



<^n+l,n+l/ 



where we have introduced the new matrix A and the 
vectors Ui and uj with: 



ui 

vi^ an 



Ui 

an 



Vi 



. (A13) 



Furt her, we need the matrices Ay defined analogously to 
dATl: 



'•'^ \ Vi-* an aij 



(A14) 



With these definitions, and with the abbreviations aij = 
det Aij and aij = det Aij , we are now able to apply Eq. 



For ttij, we make use of Eq. (A2) with n = 2, which we 
have proved above: 



\<^n+l,2 ••• <^n+l,n+ly 

(A15) 



1 



det A 



{diidij — CLiidij) ' 



(AI6) 



(A12) Inserting this result in (A15) yields a determinant with 



entries of the form anaij — anaij. We make use of 
the multi linearity of the determinant to decompose this 
expression and we obtain a sum of determinants with 
prefactors of the form aij. Eliminating zero contribu- 
tions, the resulting expression corresponds precisely to 
the Laplace-expansion of a larger determinant, and we 
finally obtain 



det Mn+i det A'' det 



^2,1 



ttl,2 
^2,2 



\<^n+l,l <^n+l,2 



<^l,n+l \ 
<^2,n+l 



<^n+l,n+l / 

(A17) 



This is the identity (A2) for n + 1. Hence we have derived 



the determinant identity for n + 1 using onl y th e identity 
for n and n = 2. By induction, the identity ( A2 ) therefore 



holds for every n G N, as it is trivial for n = 1. 
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stochastic analytic continuation. This, however is numeri- 
cally demanding and requires a very high quality of data. 
In the present case, we were unable to extract more than 
the lowest lying excitation of the dynamical charge struc- 
ture factor, which is directly connected to the charge gap. 
The higher energy spectrum showed an extremely complex 
structure which is difficult to capture with stochastic ana- 
lytic continuation. 

Strictly speaking, the reference model in the DMFT only 
has one superconducting bath, while we introduced a left 
and a right bath in the Hamiltonian ([T]). However, in the 
CTQMC, the reference model is entirely encoded in the 
bare Green's function, which can be understood as an ac- 
tion representation in the path integral formalism. The 
explicit number of the superconducting baths is therefore 
unimportant. 



